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Abstract 

A numerical method for solving Schrodinger's equation based upon a Baker-Campbell- 
Hausdorff (BCH) expansion of the time evolution operator is presented herein. The tech- 
nique manifestly preserves wavefunction norm, and it can be applied to problems in any 
number of spatial dimensions. We also identify a particular dimensionless ratio of poten- 
tial to kinetic energies as a key coupling constant. This coupling establishes characteristic 
length and time scales for a large class of low energy quantum states, and it guides the 
choice of step sizes in numerical work. Using the BCH method in conjunction with an 
imaginary time rotation, we compute low energy eigenstates for several quantum systems 
coupled to non-trivial background potentials. The approach is subsequently applied to the 
study of ID propagating wave packets and 2D bound state time development. Failures 
of classical expectations uncovered by simulations of these simple systems help develop 
quantum intuition. 

Finally, we investigate the response of a Superconducting Quantum Interference De- 
vice (SQUID) to a time dependent potential. We discuss how to engineer the potential's 
energy and time scales so that the SQUID acts as a quantum NOT gate. The notional 
simulation we present for this gate provides useful insight into the design of one candidate 
building block for a quantum computer. 
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1. Introduction 



In this article, we study the low energy behavior of a quantum system that interacts 
with an external classical background. We specifically focus upon numerically simulating 
the evolution of a system with just a few degrees of freedom. Its characteristic time scales 
are generally assumed to be short compared with the background's, yet its typical velocities 
are required to be low compared to the speed of light. The dynamics of the quantum system 
are then accurately described by the non-relativistic Schrodinger equation 



Here the "ket" denotes an abstract vector in a Hilbert space that is in one-to-one 
correspondence with the physical state of the system. 

The system's time evolution is governed by the Hamiltonian 



which is a sum of kinetic and potential energy terms. We absorb an inverse mass coefficient 
into so that the kinetic operator has unit normalization. The potential term summarizes 
the classical background's effect on the quantum system. The corresponding reaction of 
the background to the system is ignored in Schrodinger's equation. The time dependence 
of V consequently represents a fixed input rather than a dynamically determined output. 

If the potential is an arbitrary function of time, the Hamiltonian evaluated at some 
time ti typically does not commute with the Hamiltonian evaluated at another time t2- 
The formal solution to Schrodinger's equation 





H(t) =p2 + y(X,t) 




|V(t)) = U(t)|^(0)) 




then involves the unitary time evolution operator given by Dyson's formula 






^ We work with units wliere 11 = 1. We also use bold-face symbols to denote abstract quantum 
operators. 
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The T symbol in ( 1.4a ) indicates that the exponential is time ordered. As the series 
expansion in ( |1.4fc|) explicitly demonstrates, time ordering forces the temporally latest 
Hamiltonian factors to appear on the left. Given the complexity of ( p.. 41 ), it is little surprise 
that the number of time dependent problems for which U(t) can be evaluated in closed 
form is small. 

If the classical background is independent of time, the evolution operator reduces to 

U(t) = exp(-im). (1.5) 

This exponentiated sum of non-commuting and V{%) operators can be decomposed 
into a product of exponentials via the Baker-Campbell-Hausdorff (BCH) expansion: 

U(t) = e^(*^) xe"^^"^*)' [v(x),[v^(x),p2]] ^^{-itf [[v{ii),P%P^] ^-li-itf[vi^),p^]^-itv(5i)^-^tp 

(1.6) 

Since the adjoint for each exponential factor equals its matrix inverse, this decomposition 
manifestly preserves unitarity. 

Although we are ultimately interested in simulating the response of quantum systems 
to time varying backgrounds, we initially restrict our attention to problems where the 
potential is not a function of time. Expansion ( p..6| ) then provides a practical means for 
evolving state vectors on a computer. For example, one readily finds 

V'(f,t) = (x|V'(t)) = e-^*^(^)j^-i^e-**pV(^A(f ,0))j +0{t^) (1.7a) 

if just the leading exponential factors appearing on the right of ( p..6| ) are retained. Here 
JF and denote Fourier and inverse Fourier transform operations which are defined 
in Appendix A. If the next-to- leading exponentiated commutator is kept in the BCH 
decomposition of U(t), a more accurate evolution formula can be derived using operator 
identities listed in Appendix A: 

V'(x,t) = e^*'^'^(^^')-^*^(^^')^-i||e-^*p"^(7A(x,0))j +0{t^) (1.76) 

with x' = X + t^W(x). For our simulation work, we include the next-to-next-to-leading 
exponentiated nested commutator factors shown in eqn. (|1.6|) . As the resulting general- 



ization of (|1 . 76| ) is complicated and un-illuminating, we relegate it to Appendix B. 



Our numerical approach to solving Schrodinger's equation is a variant of exponentiated 
split operator methods that have been explored by several authors in the past [|],|[ . It can 
be implemented in any number of spatial dimensions via efficient Fast Fourier Transform 
(FFT) codes unlike other matrix inversion algorithms which become unwieldy beyond one 
dimension 0]. The price we pay for such generality is having to Fourier transform back 
and forth between momentum and position space. As in all split operator approaches, 
such transforming is computationally expensive. So while others have studied more com- 
plicated decompositions of U(t) that yield higher temporal order accuracy per time step 
0,^, we prefer to work with the simpler evolution formulas in (|1 . 7|) and their subleading 
generalizations which minimize the number of transforms that must be performed, i 

The outline for our article is as follows. In section 2, we identify a particular energy 
ratio as a dimensionless coupling constant that establishes characteristic length and time 
scales for several low energy quantum systems. For large values of this coupling, dimen- 
sional analysis helps guide the choice of FFT bin sizes needed to implement our BCH 
evolution algorithm. In section 3, we discuss a simple technique based upon an imagi- 
nary time rotation for projecting out low energy eigenstates from trial wavefunctions. The 
method is demonstrated to reproduce known spectra for analytically soluble ID models 
and to yield consistent eigenfunctions for analytically intractable 2D systems. In section 4, 
we investigate Gaussian wavepacket propagation on three different backgrounds. As we 
shall see, the quantum motion of a packet inside a box is, at first glance, a surprise. In 
section 5, we simulate the response of a Superconducting Quantum Interference Device 
(SQUID) to a time varying potential. By judiciously modulating SQUID parameters, this 
device can be forced to act as a quantum NOT gate. Finally we summarize our findings in 
section 6 and close with some thoughts on future simulation work for quantum computer 
design. 



2. Naive dimensional analysis 

It is generally a good practice to scale out all dimensionful parameters from 
Schrodinger's equation in order to gain qualitative insight into the quantum system it 
describes. Once dimensionful quantities have been removed, Schrodinger's equation can 

^ According to the authors of ref. |^ who surveyed several different propagation schemes, it is 
unclear whether the gain in accuracy achieved by higher order symmetric expansions of U(t) is 
offset by the considerably greater numerical effort required to implement such schemes. 
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only depend upon dimensionless ratios whose values relative to unity encode non-trivial 
physics information. In particular, these ratios act as coupling constants which govern the 
strength of the interaction between the quantum system and its classical background. 

In order to derive a dimensionless form of Schrodinger's equation, we first restore all 
dimensionful parameters within Hamiltonian ( p..2| ) and project it onto the position basis: 



We next express the position vector x as the product of some arbitrary unit of length i 
and a dimensionless vector f3 : 

The energy Kq = /2mP is subsequently scaled out from the kinetic term in (|2.1| ), and 
it is used to relate dimensionful time t to its dimensionless counterpart r: 

h 

We also rewrite the potential term as another energy Vb times a function of /3 and r: 

v{x,t) = v^u{3,T). 

The ratio 

H d'^ 

n^ — = --^ + aU{P,T) (2.2) 
with a = Vo/Kq then enters into the dimensionless Schrodinger equation 

z^V'(/3,r)=7^(r)^(/3,r). 

We work with version (|2.2| ) of the Hamiltonian in all numerical simulations. 

If the classical background is complicated, it is generally characterized by several 
length scales of different magnitudes. The potential is then implicitly a function of multiple 
length scale ratios: 

Observable expectation values can depend in complicated ways upon a and these dimen- 
sionless ratios. In this situation, naive dimensional analysis cannot be used to establish 
physical properties of low energy quantum systems such as their typical wavefunction 
spreads or oscillation frequencies. 
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But for the important special case where the background is characterized by a single 
length, dimensional considerations do fix the dependence of system scales upon the a 
coupling. For example, consider an oscillator with Hamiltonian 

n = -^ + a{Pr. (2.3) 

According to the virial theorem, the kinetic and potential energies of the oscillator's low 
lying states are balanced on average: 

{K)=p{V). (2.4) 

The states' characteristic length and time consequently scale with a as 

L~(ap)"^ T~(ap)"^, (2.5a) 

while their momentum and energy are inversely related to L and T: 

P ~ (ap) ^ E ^ (ap) ^ . (2.56) 

Note in particular that L — as a — > oo. As one would intuitively expect, an oscillator's 
wavefunction becomes more tightly concentrated about the potential's minimum as the 
potential's walls grow more steep. 

Consider next potentials such a,s U = —(3 + (3^ or U = cos (3 ^ that possess non- 
vanishing quadratic terms at their minima. In the large a limit, self-consistent dimensional 
analysis demonstrates that the scaling behavior in ( |2.5a, 6|) continues to hold with power 



p = 1. As the wavefunction becomes concentrated about the minimum point, it effectively 
senses only the leading quadratic term in the potential's expansion about the minimum. 
The low energy spectra for a large class of systems therefore look like that for a harmonic 
oscillator as a ^ oo. As an illustrative example of this general behavior, we plot in 
fig. 1 the first ten eigenvalues for periodic energy eigenstates of a one-dimensional cosine 
potential with a = 100. Though the eigenstates occur in nearly degenerate parity doublets, 
the familiar harmonic oscillator level spacing is evident in the figure. 

These elementary dimensional analysis considerations help guide the practical choice 
of FFT bin sizes needed to numerically implement the Fourier transforms in eqns. ( 1.7a, £| ) 
and their subleading generalizations: 



Low energy cosine potential spectrum 
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Fig. 1: Lowest ten eigenvalues for periodic eigenstates of a quantum system coupled to 
a U {(3) = 2 + 2 cos(2 /3) potential with strength a = 100. The eigenstates occur in 
parity even and odd pairs, but each pair's splitting is not resolvable on this plot's scale. 
The first four excited pair energies are in the ratios 2.97, 4.89, 6.76 and 8.56 relative to 
the ground state energy. Recall the corresponding harmonic oscillator energy ratios are 
3, 5, 7 and 9. 

In order to keep the position and momentum spaces on an equal footing, we set both 6(3 
and its dimensionless momentum analog dn inversely proportional to the square root of the 
total number of bins in each spatial dimension. But we modulate the bins' magnitudes 
by fractional powers of a to take into account generic systems' length and momentum 
scales. We likewise set the time step St used to evolve Schrodinger's equation proportional 
to \l\fot. These scaling choices significantly improve simulation efficiency for strongly 
coupled quantum systems. 



3. Energy eigenstate evaluation 

One of the most important properties of a quantum system is its low energy spectrum. 
It is natural to consider initially preparing a system in its ground state or in a superposition 
of a few excited states. We then want to simulate how it evolves over time and reacts to 
changes in the classical background. In order to carry out this program, we first need to 
compute the system's low energy eigenstates. 

To begin, we center C7(/3 ) about the origin /3 = 0. Translating the potential is 
especially important for simulating wavefunction evolution on periodic manifolds like 5"^ 
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(circle in ID) or x (torus in 2D). The periodicity of the FFT can be exploited 
to numerically solve Schrodinger's equation on such manifolds for states with periodic 
boundary conditions. We therefore restrict the domain for /3 to a single period in each 
spatial dimension for which U{P) is periodic. 

We next construct a trial wavefunction in position space for each low-lying energy 
state. In one spatial dimension, the parity of the ground state corresponding to a symmetric 
potential U{(3) = U{—(3) is always even. Moreover, successive non-degenerate energy states 
labeled by principle quantum number n have alternating (—1)" parity assignments. So for 
a potential which is aperiodic in /3, we take the n^^ eigenstate's trial wavefunction to be 
proportional to /S^'^^"'^'^^ exp{—/3'^). If U{/3) is periodic on the other hand, we set V'triai 
proportional to -£^Qse initial guesses are generally useful even when the 

potential does not exhibit a perfect reflection symmetry. Moreover, they correctly place 
the bulk of low energy wavefunction content in valley regions of the potential. 

Similar considerations motivate the choice of trial wavefunctions in higher spatial di- 
mensions. For example in 2D, we set the trial wavefunction for the energy eigenstate labeled 
by principle quantum numbers m and n proportional to pi^^"'^'^^ (3y^^"'^'^^U{f3j;, f3y)~^ if 
the potential is periodic in both (3j: and (3y. On the other hand, we take 'i/'triai to be a 
two-dimensional Gaussian modulated by ^(,"^™°'i2)^(nmod2) ^j^^ potential has infinite 
period. 

Once a trial wavefunction for the ground state has been selected, it can be written 
as an a priori unknown superposition of the true ground state i/^o and all other energy 
eigenstates with which it shares the same symmetry properties: 

4rL (/3 , 0) = Coi^o (/3 , 0) + ^ V, (/3 , 0) . 

j>0 

In order to distill i/'o from this infinite sum, we iteratively evolve the trial state according to 
the imaginary time evolution operator e~^^ using subleading generalizations of eqn. ( |1 . 7| ) 
with r —it: 

r) = Coe-^"-^o(/3 , 0) + J2 C,e-^^^i^,0, 0). 

i>o 

Following each timestep increment, the wavefunction is renormalized to preserve amplitude 
content. As the iteration proceeds, the trial wavefunction's excited energy state compo- 
nents become exponentially suppressed compared to the ground state term. Ultimately, 
Wriai relaxes to the true ground state 
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and its Hamiltonian expectation value approaches the genuine lowest energy eigenvalue: 



i^SUrmi^l^SUr)) Eo. 

Once the ground state is known, we select a new trial wavefunction i^H^^i, compute 
its overlap with the lowest energy state and form the difference 

^'reduced = V'iVal " ( ^0 I V'trill ) ^0 • 



We subsequently project out the first excited eigenstate from this reduced wavefunction by 
repeating the imaginary time evolution procedure described above for the ground state. By 
following this projection and evolution strategy, we can compute any energy eigenstate and 
eigenvalue so long as build-up of numerical inaccuracies does not become overwhelmingly 
large. 

As a first sanity check on this procedure, we illustrate its results in fig. 2 for a ID 
harmonic oscillator which couples with strength a = 1 to a quadratic potential. The nu- 
merically computed spectrum and eigenfunctions reproduce well-known analytic formulas. 
In order to recover dimensionful energies from the dimensionless eigenvalues shown in the 
figure, we need to multiply the latter by the energy Ko = ^h^/me'^ which we scaled out 
from Hamiltonian ( ^.2| ). Recalling Vq = ^mu'^i'^ for a simple harmonic oscillator with nat- 
ural frequency u and a = Vq/Kq, we deduce Kq = huj/2. So the n^^ state's dimensionless 
2n + 1 eigenvalue matches onto the dimensionful oscillator energy = (2n + l)hLj/2. 



Harmonic oscillator spectrum 



Initial energy eigenfunctions 
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Fig. 2: Low energy eigenvalues and eigenstates for a harmonic oscillator with coupling 
a = 1 to a quadratic potential. 




8 



As a more challenging test, we apply our eigenstate evaluation algorithm to a system 
that interacts with potential U{P) = 2 + 2 cos 2/3. In this case, the Schrodinger problem 
maps onto Mathieu's differential equation: 

+ (£;-2a) -2acos2/3 V'l/S) =0. (3.1) 



As we have have previously noted, the low energy spectrum for a ^ 1 looks like the har- 
monic oscillator's. On the other hand, known power series expressions for the eigenvalues 
of Mathieu's equation are convergent when a ^ 0{1) 0. We therefore set a = 1 and 
compare numerically derived eigenvalues with their power series counterparts in Table 1. 

Some of the discrepancies between the analytic and numerical results shown in the 
table are attributable to power series truncation. In particular, we estimate an O(10~"^) 
uncertainty in the zeroth and fourth analytic eigenvalues based upon the magnitudes of 
the last power series terms listed by Abramowitz and Stegun 0. Unfortunately, we have 
no simple way to a priori quantify errors in our numerical eigenvalues. Nevertheless, the 
agreement between the two sets of energies is overall quite good for the first 11 periodic 
eigenstates. 

As a final example, we consider a two-dimensional quantum system for which no 
analytic spectrum solution is known. It couples with unit strength to the dimensionless 
potential 

U {(3x, (3y) = 3 -|- cos 2(3y — 2 cos (3^^ cos (3y 

plotted in fig. 3. Note that U remains invariant under (3x ^ —(3x and (3y ^ —(3y refiec- 
tions, but not under fix <-> (3y exchange. The system's lowest three eigenvalues labeled by 
principle quantum numbers m and n equal Eqq = 2.59, Eqi = 3.37 and Eio = 3.79. Their 
associated eigenfunctions i/'oo, i^oi and if^io are respectively displayed in figs. 4a, 4b and 
4c. As expected, the ground state's amplitude is concentrated in valley regions of the 2D 
potential, while its phase is everywhere uniform. In contrast, the first and second excited 
states exhibit node lines along the (3y = and (3x = axes, and their amplitudes on the 
two sides of these nodal separations are 180° out of phase. Most of the excited states' 
wavefunction densities avoid the potential's mountain terrain. 

The magnitudes and relative phases of these eigenfunctions must be preserved as they 
evolve under the action of the unitary operator in ( |1 . 5| ) . We have verified that these 
stationary wavefunction conditions are indeed maintained up to small numerical errors. 
Moreover, the states' energies remain constant as time proceeds. So although we cannot 
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Table. 1: Lowest 11 energy eigenvalues for a system with periodic wavefunction bound- 
ary conditions that is coupled to potential U{f3) = 2 + 2 cos(2 f3) with strength a = 1. 
The analytic eigenvalues are based upon power series solutions to Mathieu's equation 

i- 

directly compare the results in fig. 4 with analytic expressions, our confidence in their 
validity is high. 

It is interesting to investigate time dependent interference between combinations of 
the stationary states. In fig. 5a, we display the initial superposition 

tlj{f3^,/3y) = -^(^oo(/3x,/3,) +^oi(/3x,/3y)) (3.2) 

of the ground and first excited states at r = 0. The dimensionless beat period for this 
simple combination equals Tbeat = 27r/(i?oi — Eqo) = 8.1. In fig. 5b and fig. 5c, we 
plot 'il!{(3j:, (3y) at approximately the quarter and half way points through the beat cycle. 
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ZD periodic potemrai 




Fig. 3: The two-dimensional periodic potential U {(3x^ (3y) = 3 + cos 2/3^ — 2 cos (3x cos I3y. 

Amplitude density tunneling between potential valley regions for this superposition state 
is evident in the sequence of wavefunction snapshots shown in the figure. 

One could proceed to examine the time dependence of more complicated superposition 
states. But rather than continue the bound state discussion, we turn at this point to 
investigate propagating states which exhibit qualitatively different quantum phenomena. 
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Fig. 5: (a) Initial superposition of the ground and first excited states from fig. 4. (b) 
Superposition at r ~ |;Tbeat- (c) Superposition at r ~ ^Tbeat- 
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4. Gaussian wavepacket propagation 



We start our exploration of propagating states by considering the well-known exam- 
ple of free Gaussian wavepacket evolution in one spatial dimension. The initial state is 
taken to have dimensionless mean position /3o, dimensionless standard deviation ao and 
dimensionless momentum /cq: 



V'(/9,0) 



r 1 1 


1/4 




exp 



exp [z/co/9] . 



(4.1) 



For the simple case of free propagation, Schrodinger's equation can be solved analytically 
to determine the packet's time evolution: 



V^(Ar) 



r 1 1 


1/4 


V2Txa'^{T)\ 


exp 



(/3-/3o-2fcor)- 
4a2(r) 



X (phase terms) 



where 



/ r 



(4.2) 



(4.3) 



Though the wavefunction's phase terms are exactly calculable, we choose not to explicitly 
list their ugly expressions here. 

The mean motion of the packet 



(/3(r)) = /3o + 2koT 



(4.4) 



conforms with classical intuition, i On the other hand, the packet's spreading over time 
is unexpected from particle mechanics. In order to demonstrate that the origin of this 
diffusion is not purely a classical wave effect, we restore all dimensionful parameters in 
( [1.3|) and expand the packet's width in powers of %: 



S{t) = Sq\1 1 -I- 



Am? si 



so + 0{h^). 



(4.5) 



In the classical ?i — > limit, wavepacket diffusion does not occur. So it represents a genuine 
quantum phenomenon. 

These analytic results for free Gaussian packet evolution can be reproduced by the 
numerical techniques described in the preceding sections. For example, we plot the mean 



^ Recall we have absorbed a factor of 1 /2m into our dimensionless Hamiltonian's kinetic term. 
If all dimensionful factors are restored, eqn. (|4.4D reads {x{t)) = xo +pot/m. 
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position and standard deviation of a packet with Pq = 0, ctq = l/-\/2 and fco = 4 as 
functions of time in fig. 6a. Numerical results shown in red are superposed for comparison 
on analytic results shown in green. Agreement between the two is clearly quite good. 

We consider next numerically propagating wave packets for problems where closed 
form analytic formulas are difficult to obtain. In fig. 6b, the time dependent position 
and diffusion of a packet with the same dimensionless starting point, initial width and 
momentum as that in fig. 6a are displayed. However in this second example, the packet 
encounters an infinite wall located at = 8 as represented by the dashed blue line in the 
mean position figure. It bounces off with no penetration into the energetically forbidden 
region. During the time the packet interacts with the wall, the leading and trailing portions 
of its wavefunction interfere with each other. But after the encounter is finished, the 
packet's mean motion once again looks reasonably classical, and its diffusion is essentially 
unaffected by the interaction with the wall. 

The motion of a wavepacket inside a box is much more curious. We again take the 
initial state to be the same as that in fig. 6a, and we graph its probability density in the 
first snapshot of fig. 7. Subsequent snapshots in fig. 7 depict the packet's evolution as it 
bounces between two infinite walls positioned at /? = ±8. As time proceeds, the diffusion 
of the packet becomes so pronounced that it continuously interferes with itself everywhere 
inside the box. Once the packet completely fills the energetically allowed region, it no 
longer exhibits an identifiable peak or momentum flow. Consequently, its mean position 
decays over time to zero (see fig. 6c). 

One might wonder whether this last result violates Ehrenfest's theorem which states 
that the laws of classical mechanics hold for quantum expectation values: 

™^(^) = -^(^(^))- (4-6) 

But as {V{x)) ^ V{{x)), the wavepacket's mean position does not obey Newton's law of 
motion. So while a classical particle would indefinitely ricochet off the walls of the box, the 
quantum wavepacket's mean asymptotes instead to the center. The motion of a quantum 
packet trapped inside a box is thus highly non-classical. 
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Fig. 6: Mean position and standard deviation as functions of time for a ID Gaussian 
wavepacket (a) propagating in free space, (b) bouncing off a wall and (c) bouncing inside 
a box. 
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Fig. 7: Probability density snapshots of a ID Gaussian wavepacket bouncing inside a 
box. 
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5. A notional quantum NOT gate 

We have so far investigated the interactions between simple quantum systems with 
static classical backgrounds. Now we consider modulating a potential over time and mea- 
suring the system's response. The spatial and temporal dependence of potentials corre- 
sponding to various experimental setups may be purposefully engineered to possess useful 
properties from an information theory standpoint. After a system is prepared in some ini- 
tial state, its wavefunction is intentionally manipulated via a time varying potential. The 
system's final-state wavefunction is subsequently read out and measured by the classical 
environment. This general program is currently of great interest for quantum computing 
applications. In this section, we examine the impact of a time dependent coupling upon 
one particular quantum system: the SQUID. 

Superconducting Quantum Interference Devices presently represent one of the most 
promising candidate building blocks for quantum computers. SQUID circuits can theo- 
retically support currents running clockwise and counterclockwise simultaneously. Recent 
experimental indications of such superposition states are beginning to demonstrate the va- 
lidity of quantum mechanics on macroscopic scales as well as the viability of SQUIDS 
as real- world qubits []T^-|T^. Moreover, SQUIDS are technologically attractive. These 
passive devices are small enough so that they can be mass produced on chips. Yet they 
are sufficiently large so that their quantum properties can be manipulated in a controlled 
fashion. SQUID design, fabrication and testing consequently represent active areas of 
research. 

As the quantum mechanics of SQUIDS may not be familiar to some readers, we review 
the Hamiltonian Hsquid which governs their dynamics in subsection 5.1. The potential 
that enters into Hsquid looks like a double-well within certain regions of its parameter 
space. Wavefunction localization in one of the two wells may naturally be interpreted as 
a logical "true" or "false" signal. If the couplings in Hsquid are manipulated over time, 
coherent wavefunction movement from one well to the other can be controlled. So SQUID 
quantum mechanics has enough structure for these devices to act as quantum NOT gates. 

We subsequently explore SQUID dynamics in subsection 5.2 using our simulation 
techniques. It is important to note again that the naive time dependent generalization of 
the evolution operator in (|1.5|) 

U(t)naive = exp[-i /" li{t')dt'] (5.1) 

L Jo ^ 
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differs from the true propagator in (|1.4| ) by commutator terms which enter at second order 
in the time step expansion: 

U(t)true - U(t)naive = ^ / dU / ' ^^2 [H(ti) , H(t2)] + 0(^3) . (5.2) 

These commutators render most numerical methods for solving Schrodinger's equation less 
accurate for time dependent backgrounds, and ours provides no exception to this general 
rule. But since we are more interested in rapidly gaining qualitative insight into SQUID 
experiment design than in achieving high numerical precision, we continue to use ( p..7|) and 
its subleading generalizations to propagate SQUID states interacting with time dependent 
potentials. We do, however, employ simple dynamic time step procedures to limit the 
accumulation of numerical errors. 

5.1. The SQUID Hamiltonian 

To begin the derivation of Hsquid, we sketch a cartoon of a Superconducting Quantum 
Interference Device in fig. 8. As the figure illustrates, a SQUID is simply a superconducting 
ring interrupted by a small insulator. A supercurrent of Cooper pairs fiows around the 
ring and quantum mechanically tunnels through the junction. The nonlinear relationships 
between the insulator's current and voltage were worked out by Josephson in 1962 |13[ j. 
The insulator section acts as a nonlinear circuit element, and it is referred to as a Josephson 
junction. 




Superconducting ring 



Fig. 8: An idealized SQUID pictured as a superconducting ring carrying a supercurrent 
that is interrupted by an insulator section. 
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Following Feynman [n|, we let '(/'(x ) denote the wavefunction of a low energy Cooper 
pair within the ring. In the presence of a background classical electromagnetic field, "0 
obeys a gauged Schrodinger equation. Without loss of generality, it may be decomposed 
as '0 = ^e^'^ where p = 'i/'*'^ is interpreted as Cooper pair probability density. Since i/j 
must be single-valued everywhere along the ring, its phase (/> must be a periodic function 
of X . 

Probability density p and its current counterpart J satisfy a gauged version of the 
familiar local conservation law dp/dt = —V • J. This conservation requirement implies 

j = — U* (-inw^) - ip i-ihWi)*) - 2(2e) Ai/;>1 = —\W(t)- ^ aI p 
2m I \ ml ri \ 

where A denotes the vector potential for the classical photon field. As all currents within 

a superconductor are spatially localized near its surface, J vanishes deep inside its interior. 

So 

= —A (5.3) 

/ 1/ 

everywhere along a contour which lies buried in the SQUID ring (see fig. 9). 




Fig. 9: Contour deep inside the superconducting ring along which J = 0. The total 
magnetic flux enclosed by the ring is proportional to the phase drop across the junction. 



The phase drop across the junction can be related to the total magnetic flux that passes 
through the ring's hole by integrating (|5.3|) along the contour between its endpoints: 



V(j)-ds 



2e 
2e 



2e / - 
, , A ■ ~ — (1) A ■ (is 
ri Ji n J 

^^total = (^^j^total- 



2e 



(V X A) ■ da 



(5.4) 
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The ratio $o = ^/2e represents a fundamental unit of flux. Therefore, the junction phase 
drop counts the number of "fluxons" that thread the ring. 

The total magnetic flux appearing in ( |5^ ) is typically a sum of supercurrent-generated 
internal flux plus some applied external flux: $totai = ^int + $ext- The classical expression 
for the energy associated with $int and the ring's geometrical self-inductance L motivates 
the first term in the SQUID's quantum Hamiltonian: 

W _ _ (^total - ^ext)^ _ ($o/27r)% ^2 ^ ,2 

■n-inductance 2^ 2^ 2^7 ' 



In this formula, the phase offset is related to the external flux by ( |5.4D : 0o = (27r/$o)$ext- 
A second contribution to Hsquid originates from the capacitance of the Josephson 
junction. The electric charge on the junction varies over time, for not every Cooper pair 
tunnels through the insulator. The classical energy associated with this charge motivates 
the second Hamiltonian term 



2 



TT _ Q_ _ rp 

-tlcapacitance J-^c'-^ • 

The charge operator Q = 2en effectively counts the number of pairs on the junction. When 
Hcapacitance is rewritten in terms of the dimensionless number operator n, its overall scale 
is set by the Cooper pair capacitance energy Ec = (2e)^/2C. 

The Cooper pair number and phase operators are canonically conjugate to each other, 
and they satisfy the commutation relation [n, ^] — —i. In the (j) representation, the number 
operator n —id/dcj) is seen to be the infinitesimal generator of (j) translations. Similarly, 
<j> may be interpreted as the generator of n translations. These observations are important 
for understanding the final contribution to the SQUID Hamiltonian 



■-junction 2 



Ej 

n=0 

which describes the energy needed to transport a Cooper pair across the Josephson junction 



11,111 



The quantum expression in ( |5.5| ) has no classical progenitor. Its prefactor is set by 
convention, while S+ = Xl^o 1^ ~^ ^~ ~ i^^V constitute raising and lowering 

operators. may alternatively be regarded as the finite translations e^^^ "^-^ Q±d/dn 
which alter pair number by ±1. When rewritten in terms of these exponential operators, 
the Josephson contribution reduces to its more familiar sinusoidal form 



XT. _ 

-•^junction „ 
2 



-Ej cos <f). 
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The sum of the inductance, capacitance and junction terms yields the SQUID Hamil- 
tonian 

($o/27r)2 



H 



SQUID 



2L 



The equations of motion following from this Hamiltonian 

dn i r 

^ = ^[HsQuiD,nJ 



dt 



[Hsquid, <i>] 



(5.6a) 
(5.66) 



are readily evaluated using the commutation relations listed in Appendix A. From (|5.6a| ) 
we obtain a quantum version of Kirchoff 's current conservation law 



capacitance ~l~ linductance "t~ -'■junction 



+ Ii. 







(5.7) 



with 



■•■capacitance 
linductance 



dO ^ dn 

—r- = 2e— 
dt dt 

^total $ext ($o/27r) 



/2n 



Ijunction = y—]EjS\n(t) = Ic siu ( 



(5.8a) 
(5.86) 
(5.8c) 



Conservation relation ( |5.7| ) indicates that a SQUID can be modeled by the equivalent 
circuit pictured in fig. 10. The cross appearing in the figure represents the Josephson 
junction and its supercurrent fiow ( ]5.8cD which depends upon the phase drop across the 
junction. 



1 

T 



Fig. 10: Equivalent circuit for a SQUID containing an inductor L and Josephson junction 
J. The junction's capacitance is modeled by capacitor C. 
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So long as /junction remains less than the critical value Ic = {2tv/^o)Ej, the supercur- 
rent tunnels through the insulator without dissipation. The voltage counterpart to /junction 
also depends upon the phase drop and it is obtained from the second equation of motion 
in i ^M) : 



The Josephson junction's I-V characteristics are governed by eqns. (|5.8c|) and (|5.9|) . 



We make two final comments about Hsquid- Firstly, it is customary to introduce the 
dimensionless parameter 

which counts the number of "fluxons" that thread the SQUID ring when it carries critical 
current Ic- When parameter L is eliminated in favor of Pl, the Hamiltonian in the (p 
representation reduces to 

HsQvm = -Ea^ + E, [ ^^:/°^' - cos . (5.10) 

Secondly, we scale out capacitance energy Ec from //^squid for simulation purposes. The 
dimensionless SQUID Hamiltonian then takes the same form as the template in (|2.2|) 



n = -—-r + aU{<P) (5.11) 



with coupling constant a = Ej/Ec and potential 

U = — — cos(j). (5.12) 

5.2. SQUID response to a time dependent background^ 

The SQUID Hamiltonian is a function of the two free parameters and 4>o in addition 
to the overall coupling a = Ej/Eq. If naive dimensional analysis intuition is to hold, (3^ 
and (f)Q must both be of order unity. Their precise values control the potential's shape 
and thereby strongly influence one's ability to extract useful information from low lying 



^ After work on this article was completed, we learned that a similar analysis of SQUID 



response to a time dependent pulse was recently reported in ref. [17|. Where overlap exists, there 



is generally good agreement between the findings of ref. [17| and the results independently derived 
in this subsection. 
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n I 1 1 1 1 1 1 

-3%I2 -V. -3t/2 a 1112 V. 3%/2 

Fig. 11: Dimensionless SQUID potential corresponding to parameters (3l = (j^o = and 
coupling a = Ej/Eq = 10. The colored horizontal lines denote the system's lowest four 
energy eigenvalues. The splitting between the ground and first excited state energies is 
not resolvable on this plot's energy scale. 

quantum states. We choose to G.x (3l = 4>o = so that the SQUID potential looks like a 
double- well centered about (p = (po- For simulation as well as visualization purposes, we 
also initially set a = 10. U{(fi) is plotted for these parameter choices in fig. 11. 

The Hamiltonian remains invariant under a discrete (p — n —{(p — n) symmetry, and 
the system's energy eigenstates are even and odd with respect to this parity operation. As 
the horizontal energy lines in fig. 11 demonstrate, the low lying parity partners are nearly 
degenerate: 

E^^^ = 13.8916 E^r^ = 13.8960 

(+) (-) ^^-^^^ 

El^' = 17.7426 E[ ' = 17.9174. 

These energy eigenvalues establish relevant time scales for the design of a quantum NOT 
gate. For example, suppose the SQUID is initially prepared in the state ^p = {ipQ~^^ + 
ipQ ^)/y/2 which is localized in the positive well and corresponds to a counterclockwise 
current. In half a beat period 0.5Tbeat = ^/{Eq ^ — Eq~^^) = 714, the system naturally 
tunnels into the negative well, and the current flows in a clockwise direction. Intentional 
switching of the SQUID's wavefunction must obviously be performed and measured on a 
much shorter time scale. H 



In this notional NOT gate example, we make no attempt to model decoherence effects which 

24 



o 

II 



II 



30 



I 1 1 


1 t ■ ill 
1 I I If 


V 

1 


J 1 .b .1 . 




Fig. 12a 
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Fig, 12c 



are present in any real- world quantum experiment. But it is important to realize that the deco- 
herence time could well set a much more stringent upper bound on when a measurement must be 
performed than does O.STbeat. 
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II 




rifts 



Fig. 12d 



Fig. 12: Snapshots of the dimensionless potential for a SQUID with (3l = (t>o = and 
variable a. The light blue horizontal lines represent the system's time dependent energy. 
The evolution of the SQUID's probability density illustrates the system's response to 
the background. 

The energy of the initial state as well as its probability density are graphed in fig. 12a. 
As can be seen in the figure, the horizontal light blue line representing the state's energy 
lies significantly below the local center maximum of the SQUID potential. Probability 
density oscillation between the wells is classically forbidden, and it only slowly proceeds 
quantum mechanically. In order to increase the frequency of wavefunction beating, we 
need to lower the potential barrier. This reduction may be implemented by diminishing 
the coupling constant. 

Modulating a SQUID's a = Ej/Ec value can be experimentally achieved by subjecting 
the Josephson junction to an external magnetic field. The rate at which a is modified is 
constrained by the information theory requirement that high energy SQUID states not be 
unduly excited. In particular, we do not want the potential's time dependence to introduce 
excitations which are much more energetic than the separation 



e[+'^ + E\ 



AE 



10 



(-) 



-^0 + -^0 



between the zeroth and first pairs of SQUID eigenstates. This second relevant energy scale 
provides an order-of-magnitude estimate rtrans ~ 27r/Ai?io = 1.6 for the time during which 
a should transition from its large initial value to a smaller intermediate size. 
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We next use our simulator to determine an intermediate coupling for which well- 
hopping becomes classically allowed. We find that if a is reduced from 10 to 0.4 according 
to the schedule shown in fig. 13, the resulting SQUID energy lies slightly above the modified 
potential's barrier. The dashed vertical red lines in the figure depict the transition time 
Ttrans- The potential, system energy and system probability density after the transition 
is complete are illustrated in fig. 12b. As we see from the eigenenergies for the a = 0.4 
potential 



= 0.7377 e^~^ 



0.9834 



(5.14) 



- 1.5514 e\ ' = 2.1090, 

the splitting between the ground and first excited state is two orders of magnitude larger 
than that for its a = 10 counterpart. So the SQUID's probability density migrates nearly 
100 times faster from one well to the other following the transition. 



Coupling constant time dependence 



S 6 




Fig. 13: SQUID coupling time dependence which induces NOT gate behavior: 



a(r) = 10 + 



10 - 0.4 



r/ r-h\ /r- 13.8 

err ( — \ — err 



0.8^2^ 



1.6^2 



After an interval approximately equal to one half of the new beat period ^-^T^^eat ~ 
7r/(£Q ■* — Eq^'') = 12.8, most of the SQUID's probability density resides within the negative 
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well (see fig. 12c). We then want to raise the barrier back to its original height in order to 
trap the bulk of the wavefunction in its new location. The separation 

^^10 — 7^ 7^ 



between the zeroth and first pairs of a = 0.4 SQUID eigenstates again provides an order- 
of-magnitude upper bound for a reasonable second transition time r^'j-^ns ~ Syr/Aeio = 6.5. 
After running several simulation tests, we find that taking r/j.a,ns = 3.2 leads to successful 
wavefunction capture with a final SQUID energy -Egnai = 14.3 that is quite close to its 
original value -Einit = 13.9. The timing and duration of the potential's return to its original 
form are depicted by the green lines in fig. 13. The system's probability density after the 
entire coherent operation is over is shown in fig. 12d. 

A measurement made of the SQUID 's final current sense is highly likely to be opposite 
to that of its initial direction. So driving the a = Ej/Ec coupling according to a schedule 
like that in fig. 13 represents one possible way to implement a quantum NOT gate. It is 
interesting to note that all quantum computing logic can be built up from NOT gates as 
well as CNOT (controlled not) gates. SQUID implementations of CNOT have also been 
considered in the recent literature [ p!8| . 



6. Conclusion 

In this article, we have developed a new approach for numerically solving Schrodinger's 
equation, and we have used it to analyze the low energy behavior of several quantum 
systems. Our numerical algorithms are based upon a Baker-Campbell-Hausdorff expansion 
of the time evolution operator which manifestly preserves unitarity and works in any 
number of spatial dimensions. We have also identified a ratio of characteristic potential 
to kinetic energies as a key coupling constant a that fixes the strength of the interaction 
between a quantum system and its classical background. For problems where a represents 
the only dimensionless parameter whose value significantly exceeds unity, dimensional 
analysis establishes relevant length and time scales for low energy states. It is important to 
take the a dependence of these physical scales into account when numerically integrating 
Schrodinger's equation. 

We have applied our numerical techniques to compute energy eigenvalues and eigen- 
states for a number of examples. The methodology has been validated by reproducing 
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known answers for analytically soluble problems and computing consistent results for an- 
alytically intractable models. We have also simulated the evolution of bound states and 
propagating packets. Various quantum phenomena such as potential barrier tunneling, 
stationary state interference, and wavefunction diffusion have been graphically displayed. 
These illustrative examples help develop quantum intuition, especially when they clash 
with classical expectations. 

Finally, we have demonstrated that our approach can be usefully applied to systems 
interacting with time dependent potentials. In particular, it provides a valuable tool for 
designing quantum information devices. Theoretical simulation of quantum circuits can 
help guide ongoing and future experiments which are important from both fundamental 
physics and technological standpoints. Indeed, this last point provided the impetus for 
this work. 

In closing, we briefly mention a few applications of the methods and insights presented 
in this article which we believe will be interesting to pursue. One current outstanding chal- 
lenge in quantum computing is to build a robust measuring apparatus that can read out a 
qubit's state without prematurely destroying its quantum information. SQUID technology 
is being pursued for not only qubit construction but also measurement operations. Pre- 
liminary simulations of such SQUID measuring circuits evince many of the same quantum 
issues as the more elementary examples discussed in this paper: a = Ej/Eq values which 
are orders of magnitude larger than unity, initial wavefunctions which are concentrated in 
deep potential wells, temporal pulse shaping which is needed to avoid undue excitations, 
time scales for wavepacket mean motion which are longer than those for packet spreading, 
etc. Theoretical simulation and experimental implementation of SQUID measurement cir- 
cuits will be reported elsewhere [0. But clearly, many of the ideas discussed here can 
be usefully applied to exciting problems that lie at the quantum information processing 
frontier. 
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Appendix A. Position and momentum space conventions 

We collect together in this appendix various useful position and momentum eigenstate 
and operator identities in n spatial dimensions. We follow Feynman's approach to keeping 
track of where factors of 2n appear in these relations. One factor of 2n accompanies every 
momentum space delta function in each spatial dimension, while l/2n accompanies every 
momentum space measure factor. No other sources of 27r enter into the identities below. 

Eigenstate orthonormality relations: 



{x\x') 
ip\p') 



Eigenstate coui|)lcteuess relations: 




Delta function identities: 




Position and momentum eigenstate overlap: 





Fourier transform conventions: 




Here i(^{x) = {x\i(^) denotes a position space wavefunction, while i(^(p) = {p\i(^) repre- 
sents its momentum space counterpart. 
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Commutation relations: I 



[P,X] =-i 
[P,y(X)] =-zVy(X) 

[P^ ViX)] = 2{-iVV{±)] ■ P + (-zV)V(X) 

Appendix B. Explicit but approximate O(t^) wavefunction evolution formula 

After expanding all commutators in the time evolution operator U(t) shown in ( |1.6|) , 
we let it act upon an initial state vector |i/'(0)): 

For computational speed purposes, we want to derive a wavefunction evolution formula 
that involves just one Fourier transform and one inverse transform. So we approximate 
E^.P.V^yP, ^ CP 2 and Ej>fc(5fc5j^)PjPfc ~ E,>fcC',fcP,Pfc for some potential- 
dependent constants C and Cjk- We then project the state vector onto the position basis 
to deduce 

where x ' = x + t^VV{x). 



^ Recall that the commutator of two hermitian operators always equals an anti-hermitian 
combination of operators. 
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